function [T, Q, pcount,zcount] = AR_Q(Ts, ana_sig, P_min, P_max, q_min, q_max,Ttarget,Qtarget)
%AR or Sompi method
%josh crozier 2020
%adapted from Lesage 2008

%Minimum and maximum number of poles of the estimated ARMA filters. 
%The minimum number of poles should be approximately twice the number of 
%clear spectral peaks.
% P_min = 4;
% P_max = 8;
%Minimum and maximum number of zeros of the estimated ARMA filters. 
%A small number of zeroes is useful to modelize noisy signals. 
%If q is null only AR filters are calculated.
% q_min = 0;
% q_max = 4;

% Autocorrelation of the selected signal
autocor=xcorr(ana_sig);

% Normalization of the autocorrelation
autocor=autocor/max(autocor);

freal = zeros((q_max-q_min+1)*(P_max-P_min+1),1);
g = zeros((q_max-q_min+1)*(P_max-P_min+1),1);
Q = zeros((q_max-q_min+1)*(P_max-P_min+1),1);
pcount = NaN(size(freal));
zcount = NaN(size(freal));
n_pole = 1;

for q=q_min:q_max
  for p=P_min:P_max
   
    filter_ar = get_ar(autocor,p,q);
%     filter_poly_ar = poly(filter_ar);

    % Calculates the corresponding complex frequencies
    
    x=real(filter_ar);
    y=imag(filter_ar);
    r=sqrt(x.^2+y.^2);
    theta=atan2(y,x);

    h=find(theta<=0 | r>1.00);
    r(h,:)=[];
    %theta(theta<=0)=[];
    theta(h)=[];
    nb_p = length(theta);
    % Frequency and quality factor
    f = theta./(Ts.*2*pi);
    qf = theta./(2*(1-r));
    freal(n_pole:n_pole+nb_p-1) = f;
    pcount(n_pole:n_pole+nb_p-1) = p;
    zcount(n_pole:n_pole+nb_p-1) = q;
    g(n_pole:n_pole+nb_p-1)=-(1-r)./(Ts*2*pi);
    
    Q(n_pole:n_pole+nb_p-1)=qf;
%     Q(n_pole:n_pole+nb_p-1) = -2*pi*freal(n_pole:n_pole+nb_p-1)./(g(n_pole:n_pole+nb_p-1));

    n_pole = n_pole + nb_p;
    % Write the f and Q of all the poles
%     str = ['p = ' num2str(p,'%3.0f') '  q = ' num2str(q,'%3.0f')];
%     fprintf(fid_fq,'%s\n',str);
%     for ip = 1:nb_p
%         str = ['f = ' num2str(f(ip),'%5.2f') ' Hz   Q = ' num2str(qf(ip),'%5.2f')];
%         fprintf(fid_fq,'%s\n',str);
%         if qf(ip) < 0
%             msg = ['p = ' num2str(p,'%3.0f') '  q = ' num2str(q,'%3.0f')...
%                 ' r = ' num2str(r(ip),'%5.2f') ' Q = ' num2str(qf(ip),'%5.2f')]
%         end
%     end
  end
end

%Graphics of the complex frequency plane
% f_g_graphic(Ts, freal, g)
% f_Q_graphic(Ts, freal, Q, Ttarget,Qtarget)
% time_graphic(hObject, eventdata, handles);
T = 1./freal;

end %end function

